Note

Please, click on this link to get back to the Home Page.

Here I will compare results of packages HMSC, gllvm, ecoCopula and EMtree.

Packages & data

# Packages

## Models

library(Hmsc)
library(ecoCopula)
library(PLNmodels)
library(EMtree)
library(gllvm)

## Graph
library(ggraph)
library(tidygraph)

## Plots
library(lattice)
library(corrplot)
library(gclus)
library(viridis)

# Data

data("spider")

sp_names <- colnames(spider$abund)

abund <- as.matrix(spider$abund)

env_sc <- data.frame(scale(spider$x))

## seed

set.seed(42)

## Number of chains for MCMC

nChains <- 3

## Number of cores for parallel computing

if (parallel::detectCores() > 2 & parallel::detectCores() > nChains) {
    
    nParallel <- nChains
    
} else {
    
    nParallel <- 1
    
}

HMSC

First, let’s try HMSC. We don’t have informations about sampling design, let’s assume that each site have been visited once.

studyDesign <- data.frame(sample = as.factor(1:nrow(abund)))

rL <- HmscRandomLevel(units = studyDesign$sample)

Now, we will assume that there is no interactions between covariables.

model_hmsc <- Hmsc(Y = abund, XFormula = ~soil.dry + bare.sand + fallen.leaves + 
    moss + herb.layer + reflection, XData = env_sc, studyDesign = studyDesign, ranLevels = list(sample = rL))

Then we set up the MCMC chains.

if (.Platform$OS.type != "windows") {
    thin <- 10
    samples <- 500
    transient <- 5 * thin
    verbose <- 1
} else {
    thin <- 10
    samples <- 5500
    transient <- 50 * thin
    verbose <- 1
}

And run the model.

model_hmsc <- sampleMcmc(model_hmsc, thin = thin, samples = samples, transient = transient, 
    nChains = nChains, verbose = 1, nParallel = nParallel)

Then we diagnostic the model.

model_hmsc_post <- convertToCodaObject(model_hmsc)

par(mfrow = c(2, 2))

hist(effectiveSize(model_hmsc_post$Beta), main = "ess(beta)")
hist(gelman.diag(model_hmsc_post$Beta, multivariate = FALSE)$psrf, main = "psrf(beta)")
hist(effectiveSize(model_hmsc_post$Omega[[1]]), main = "ess(omega)")
hist(gelman.diag(model_hmsc_post$Omega[[1]], multivariate = FALSE)$psrf, main = "psrf(omega)")

par(mfrow = c(1, 1))

coda::traceplot(model_hmsc_post$Beta)

Not so great for chain convergence…

Let’s see the correlation with the environnemtal variables.

postBeta = getPostEstimate(model_hmsc, parName = "Beta")
plotBeta(model_hmsc, post = postBeta, param = "Support", supportLevel = 0.95)

There is no correlation with the species and the bare sediment and reflexion.

OmegaCor <- computeAssociations(model_hmsc)
supportLevel <- 0.95

toPlot <- ((OmegaCor[[1]]$support > supportLevel) + (OmegaCor[[1]]$support < (1 - 
    supportLevel)) > 0) * OmegaCor[[1]]$mean

corrplot(toPlot, method = "color", col = colorRampPalette(c("blue", "white", "red"))(200), 
    title = paste("random effect level:", model_hmsc$rLNames[1]), mar = c(0, 0, 1, 
        0))

So there is no inter-specific interactions.

gllvm

First, let’s choose the number of latent variables.

criterias <- NULL

for (i in 0:5) {
    
    fiti <- gllvm(abund, env_sc, family = "negative.binomial", num.lv = i, formula = abund ~ 
        soil.dry + bare.sand + fallen.leaves + moss + herb.layer + reflection)
    
    criterias[i + 1] <- summary(fiti)$AICc
    
    names(criterias)[i + 1] = i
}

print(criterias)

       0        1        2        3        4        5 
1137.907 1131.179 1170.915 1200.526 1229.728 1258.644 

Based on \(AIC\), We will take 1 latent variable. Now we are fitting the model.

model_gllvm <- gllvm(abund, env_sc, family = "negative.binomial", num.lv = 1, formula = abund ~ 
    soil.dry + bare.sand + fallen.leaves + moss + herb.layer + reflection)

Let’s see some diagnostic plots.

par(mfrow = c(3, 2), mar = c(4, 4, 2, 1))
plot(model_gllvm, var.colors = 1)

Effect of environnemtal variables.

coefplot(model_gllvm, cex.ylab = 0.7, mfrow = c(1, 1))

Species interactions

cr0 <- getResidualCor(model_gllvm)

corrplot(cr0[order.single(cr0), order.single(cr0)], diag = FALSE, type = "lower", 
    method = "square", tl.cex = 0.8, tl.srt = 45, tl.col = "red")

Besides one specie (Arctperi), all have negative residual correlation.

ecoCopula

GCGM_mod <- manyglm(mvabund(abund) ~ as.matrix(env_sc))
GCGM_graph <- cgr(GCGM_mod)

plot(GCGM_graph, pad = 1)

It’s making a small graph, some species do not seems to interact.

EMTree

We’ll start with the same model

EMT_model <- PLN(abund ~ as.matrix(env_sc))

 Initialization...
 Adjusting a PLN model with full covariance model
 Post-treatments...
 DONE!

p <- ncol(abund)  # Number of nodes

EMT_model_res <- EMtree(EMT_model, maxIter = 10, plot = TRUE)

Likelihoods: 17.61764 , 17.6573 , 17.65734 , 

Convergence took 0.27 secs  and  3  iterations.
Likelihood difference = 4.241987e-05 
Betas difference = 8.172318e-08

str(EMT_model_res)

List of 5
 $ edges_prob  : num [1:12, 1:12] 0 0.2777 0.1111 0.1336 0.0908 ...
 $ edges_weight: num [1:12, 1:12] 0 0.00761 0.00756 0.00757 0.00755 ...
 $ logpY       : num [1:3] 17.6 17.7 17.7
 $ maxIter     : num 3
 $ timeEM      : 'difftime' num 0.26788592338562
  ..- attr(*, "units")= chr "secs"

edges_prob <- EMT_model_res$edges_prob

edges_prob[edges_prob < 2/p] <- 0  # Remove links too improbable

EMtree::draw_network(edges_prob, nodes_label = sp_names, pal = "dodgerblue3", layout = "nicely", 
    curv = 0.1)

$G

$graph_data
# A tbl_graph: 12 nodes and 20 edges
#
# An undirected simple graph with 2 components
#
# Node Data: 12 x 8 (active)
    btw bool_btw bool_deg   deg title name     label      finalcolor
  <dbl> <lgl>    <lgl>    <dbl> <chr> <chr>    <chr>      <lgl>     
1     8 FALSE    TRUE         4 ""    Alopacce "Alopacce" FALSE     
2     0 FALSE    TRUE         4 ""    Alopcune "Alopcune" FALSE     
3     0 FALSE    TRUE         2 ""    Alopfabr "Alopfabr" FALSE     
4     0 FALSE    TRUE         2 ""    Arctlute "Arctlute" FALSE     
5     0 FALSE    FALSE        0 ""    Arctperi ""         FALSE     
6     0 FALSE    TRUE         2 ""    Auloalbi "Auloalbi" FALSE     
# … with 6 more rows
#
# Edge Data: 20 x 6
   from    to weight btw.weights neibs title
  <int> <int>  <dbl>       <dbl> <lgl> <chr>
1     1     2  0.278        1.53 FALSE ""   
2     1     8  0.369        1.31 FALSE ""   
3     1     9  0.362        1.32 TRUE  ""   
# … with 17 more rows